﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// The class represents a mathematic matrix for non complex numbers.
    /// </summary>
    [Serializable]
    public class Matrix : AbstractMatrix, IMatrix
    {
        /// <summary>
        /// The precision of the results for solving a linear system.
        /// </summary>
        private const double Precision = 0.00000001;

        /// <summary>
        /// The number of rows and columns of the matrix.
        /// </summary>
        private readonly int rows;

        /// <summary>
        /// The number of rows and columns of the matrix.
        /// </summary>
        private readonly int columns;

        /// <summary>
        /// The coefficients of the matrix.
        /// </summary>
        private readonly double[,] matrixData;

        /// <summary>
        /// Initializes a new instance of the <see cref="Matrix"/> class.
        /// </summary>
        /// <param name="rows">The number of matrix rows.</param>
        /// <param name="columns">The number of matrix columns.</param>
        public Matrix(int rows, int columns)
        {
            if (rows < 0)
            {
                throw new ArgumentException("rows < 0", "rows");
            }

            if (columns < 0)
            {
                throw new ArgumentException("columns < 0", "columns");
            }

            this.columns = columns;
            this.rows = rows;
            this.matrixData = new Double[rows,columns];
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Matrix"/> class.
        /// </summary>
        /// <param name="matrix">The matrix as a two dimensional array.</param>
        public Matrix(Double[,] matrix)
        {
            if (matrix == (double[,]) ((object) null))
            {
                throw new ArgumentNullException("matrix");
            }

            this.rows = matrix.GetLength(0);
            this.columns = matrix.GetLength(1);
            this.matrixData = new double[matrix.GetLength(0),matrix.GetLength(1)];

            for (int i = 0; i < matrix.GetLength(0); i++)
            {
                for (int j = 0; j < matrix.GetLength(1); j++)
                {
                    this.matrixData[i, j] = matrix[i, j];
                }
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Matrix"/> class.
        /// </summary>
        /// <param name="matrix">The matrix to clone.</param>
        public Matrix(Matrix matrix)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            this.rows = matrix.Rows;
            this.columns = matrix.Columns;
            this.matrixData = (double[,]) matrix.MatrixData.Clone();
        }

        /// <summary>
        /// Gets the number of matrix columns.
        /// </summary>
        public int Columns
        {
            get { return columns; }
        }

        /// <summary>
        /// Gets the number of matrix rows.
        /// </summary>
        public int Rows
        {
            get { return rows; }
        }

        /// <summary>
        /// Gets the matrix data in form of an array.
        /// </summary>
        public double[,] MatrixData
        {
            get { return matrixData; }
        }

        /// <summary>
        /// Creates the unit matrix of the current matrix. The new matrix will overwrite the
        /// old one.
        /// </summary>
        public void CreateIdentityMatrix()
        {
            for (int i = 0; i < this.rows; i++)
            {
                if (ExtendedMath.Fabs(this.matrixData[i, i]) > Matrix.Precision)
                {
                    for (int k = i + 1; k < this.columns; k++)
                    {
                        this.matrixData[i, k] /= this.matrixData[i, i];
                    }

                    this.matrixData[i, i] = 1;
                }
            }
        }

        /// <summary>
        /// Solves the linear system of the matrix by using the Gauss elimination method. 
        /// To solve the system of equation, the matrix has to be am extended coefficient matrix.
        /// The new matrix will overwrite the old one.
        /// </summary>
        /// <returns>The results of the system solve.</returns>
        public Double[] SolveLinearSystem()
        {
            this.CreateUpperTriangularMatrix();
            this.CreateDiagonalMatrix();
            this.CreateIdentityMatrix();

            Double[] result = new Double[this.rows];

            for (int i = 0; i < this.rows; i++)
            {
                result[i] = this.matrixData[i, this.columns - 1];
            }

            return result;
        }

        /// <summary>
        /// Creates the upper triangular matrix. The new matrix will overwrite the
        /// old one.
        /// </summary>
        public void CreateUpperTriangularMatrix()
        {
            int reminder;
            double tempuri;

            for (int i = 0; i < this.rows - 1; i++)
            {
                reminder = this.GetRowWithBiggestElement(i);

                if (reminder != i)
                {
                    this.SwapRow(i, reminder);
                }

                if (ExtendedMath.Fabs(this.matrixData[i, i]) > Matrix.Precision)
                {
                    for (int k = i + 1; k < this.rows; k++)
                    {
                        if (ExtendedMath.Fabs(this.matrixData[k, i]) > Matrix.Precision)
                        {
                            tempuri = this.matrixData[k, i] / this.matrixData[i, i];

                            for (int j = i + 1; j < this.columns; j++)
                            {
                                this.matrixData[k, j] -= tempuri * this.matrixData[i, j];
                            }
                        }

                        this.matrixData[k, i] = 0.0;
                    }
                }
                else
                {
                    this.matrixData[i, i] = 0.0;
                }
            }
        }

        /// <summary>
        /// Creates the diagonal matrix. The new matrix will overwrite the 
        /// old one.
        /// </summary>
        public void CreateDiagonalMatrix()
        {
            double tempuri;

            for (int i = this.rows - 1; i >= 1; i--)
            {
                if (ExtendedMath.Fabs(this.matrixData[i, i]) > Matrix.Precision)
                {
                    for (int k = i - 1; k >= 0; k--)
                    {
                        if (ExtendedMath.Fabs(this.matrixData[k, i]) > Matrix.Precision)
                        {
                            tempuri = this.matrixData[k, i] / this.matrixData[i, i];

                            for (int j = i + 1; j < this.columns; j++)
                            {
                                this.matrixData[k, j] -= tempuri * this.matrixData[i, j];
                            }
                        }

                        this.matrixData[k, i] = 0;
                    }
                }
            }
        }

        /// <summary>
        /// Calculates the main diagonal product of the current matrix.
        /// </summary>
        /// <returns>The main diagonal product of the current matrix.</returns>
        public double CalculateMainDiagonalProduct()
        {
            if (this.IsSquareMatrix)
            {
                double product = 1;

                for (int i = 0; i < this.rows; i++)
                {
                    product *= this.matrixData[i, i];
                }

                return product;
            }

            throw new IllegalArithmeticException("Cannot calculate the main diagonal product of a non square matrix.");
        }

        /// <summary>
        /// Calculates the determinant of the current matrix.
        /// </summary>
        /// <returns>The determinant of the current matrix.</returns>
        public double CalculateDeterminant()
        {
            Matrix tempuri = new Matrix(this.matrixData);

            tempuri.CreateUpperTriangularMatrix();

            return -tempuri.CalculateMainDiagonalProduct();
        }

        /// <summary>
        /// Creates the sign inverted matrix of the current instance.
        /// </summary>
        /// <returns>The sign inverted matrix of the current instance.</returns>
        public Matrix CreateSignInvertedMatrix()
        {
            Matrix result = new Matrix(this.rows, this.columns);
            ;

            for (int i = 0; i < this.rows; i++)
            {
                for (int j = 0; j < this.columns; j++)
                {
                    result.SetValueAtPosition(i, j, -this.matrixData[i, j]);
                }
            }

            return result;
        }

        /// <summary>
        /// Returns the row with biggest element.
        /// </summary>
        /// <param name="row">The row.</param>
        /// <returns>The index of the biggest element in a matrix row.</returns>
        private int GetRowWithBiggestElement(int row)
        {
            int result = row;
            double tempuri = 0;

            tempuri = ExtendedMath.Fabs(this.matrixData[row, row]);

            for (int i = row + 1; i < this.rows; i++)
            {
                if (ExtendedMath.Fabs(this.matrixData[i, row]) > tempuri)
                {
                    tempuri = ExtendedMath.Fabs(this.matrixData[i, row]);
                    result = i;
                }
            }

            return result;
        }

        /// <summary>
        /// Transposes the current matrix.
        /// </summary>
        /// <returns>The transposed matrix.</returns>
        public Matrix Transpose()
        {
            Matrix result = new Matrix(this.columns, this.rows);

            for (int i = 0; i < this.rows; i++)
            {
                for (int j = 0; j < this.columns; j++)
                {
                    result.SetValueAtPosition(new MatrixPosition(j, i),
                                              this.GetValueAtPosition(new MatrixPosition(i, j)));
                }
            }

            return result;
        }

        /// <summary>
        /// Swaps the row.
        /// </summary>
        /// <param name="row1">The first row.</param>
        /// <param name="row2">The row which should be swap with the first row.</param>
        private void SwapRow(int row1, int row2)
        {
            double tempuri = 0;

            for (int i = 0; i < this.columns; i++)
            {
                tempuri = this.matrixData[row1, i];
                this.matrixData[row1, i] = this.matrixData[row2, i];
                this.matrixData[row2, i] = tempuri;
            }
        }

        /// <summary>
        /// Gets the value at the specified position of the matrix.
        /// </summary>
        /// <param name="position">The specified position where the value should be set.</param>
        /// <returns>The value at the specified position.</returns>
        public Double GetValueAtPosition(MatrixPosition position)
        {
            if ((position.Row < this.matrixData.GetLength(0)) && (position.Column < this.matrixData.GetLength(1)))
            {
                return this.matrixData[position.Row, position.Column];
            }

            throw new IllegalPositionException("Matrix position does not exist");
        }

        /// <summary>
        /// Gets the value at the specified position of the matrix.
        /// </summary>
        /// <param name="row">The row index where the value should be set.</param>
        /// <param name="column">The column index where the value should be set.</param>
        /// <returns>The value at the specified position.</returns>
        public Double GetValueAtPosition(int row, int column)
        {
            if ((row < this.matrixData.GetLength(0)) && (column < this.matrixData.GetLength(1)))
            {
                return this.matrixData[row, column];
            }

            throw new IllegalPositionException("Matrix position does not exist");
        }

        /// <summary>
        /// Sets a specified value at a specified position of the matrix.
        /// </summary>
        /// <param name="position">The position where the value should be set.</param>
        /// <param name="value">The number of the element.</param>
        public void SetValueAtPosition(MatrixPosition position, Double value)
        {
            if ((position.Row < this.matrixData.GetLength(0)) && (position.Column < this.matrixData.GetLength(1)))
            {
                this.matrixData[position.Row, position.Column] = value;
            }
            else
            {
                throw new IllegalPositionException("Matrix position does not exist");
            }
        }

        /// <summary>
        /// Sets a specified value at a specified position of the matrix.
        /// </summary>
        /// <param name="row">The row index of the new value.</param>
        /// <param name="column">The column index of the new value.</param>
        /// <param name="value">The number of the element.</param>
        public void SetValueAtPosition(int row, int column, Double value)
        {
            if ((row < this.matrixData.GetLength(0)) && (column < this.matrixData.GetLength(1)))
            {
                this.matrixData[row, column] = value;
            }
            else
            {
                throw new IllegalPositionException("Matrix position does not exist");
            }
        }

        /// <summary>
        /// Implements the operator *.
        /// </summary>
        /// <param name="matrix">The matrix.</param>
        /// <param name="vector">The vector for the multiplication.</param>
        /// <returns>The result of the operator.</returns>
        public static GeneralVector operator *(Matrix matrix, GeneralVector vector)
        {
            if (vector == (GeneralVector) null)
            {
                throw new ArgumentNullException("vector");
            }

            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (vector.Count != matrix.Columns)
            {
                throw new ArithmeticException(
                    "The number of vector rows and the number of matrix columns have to be even.");
            }

            GeneralVector result = new GeneralVector(matrix.Rows);

            for (int i = 0; i < matrix.Rows; i++)
            {
                double tempuri = 0;

                for (int j = 0; j < matrix.Columns; j++)
                {
                    tempuri += matrix.GetValueAtPosition(i, j) * vector[j];
                }

                result[i] = tempuri;
            }

            return result;
        }

        /// <summary>
        /// Implements the operator *.
        /// </summary>
        /// <param name="matrix">The matrix.</param>
        /// <param name="number">The scalar for the multiplication.</param>
        /// <returns>The result of the operator.</returns>
        public static Matrix operator *(Matrix matrix, Double number)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            Matrix result = matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    result.SetValueAtPosition(new MatrixPosition(i, j),
                                              result.GetValueAtPosition(new MatrixPosition(i, j)) * number);
                }
            }

            return result;
        }

        /// <summary>
        /// Appends a identity matrix to the specified matrix.
        /// </summary>
        /// <param name="a">The source matrix.</param>
        /// <returns>The source matrix with the append identity matrix.</returns>
        public static Matrix AppendIdentityMatrix(Matrix a)
        {
            if (a == (Matrix) null)
            {
                throw new ArgumentNullException("a");
            }

            Matrix result = new Matrix(a.Rows, a.Columns + a.Rows);

            for (int i = 0; i < a.Rows; i++)
            {
                for (int j = 0; j < a.Columns; j++)
                {
                    result.SetValueAtPosition(new MatrixPosition(i, j),
                                              a.GetValueAtPosition(new MatrixPosition(i, j)));
                }
            }

            for (int i = 0; i < a.Rows; i++)
            {
                result.SetValueAtPosition(new MatrixPosition(i, i + a.Columns), 1);
            }

            return result;
        }

        /// <summary>
        /// Unions the specified matrices. The two matrices have to have the same number of
        /// rows.
        /// </summary>
        /// <param name="a">The first matrix.</param>
        /// <param name="b">The second matrix.</param>
        /// <returns>The union matrix.</returns>
        public static Matrix UnionMatricesVertical(Matrix a, Matrix b)
        {
            if (a.Rows != b.Rows)
            {
                throw new ArgumentException("The number of rows in both matrices has to be even.");
            }

            Matrix result = new Matrix(a.Rows, a.Columns + b.Columns);

            for (int i = 0; i < a.Rows; i++)
            {
                for (int j = 0; j < a.Columns; j++)
                {
                    result.SetValueAtPosition(new MatrixPosition(i, j),
                                              a.GetValueAtPosition(new MatrixPosition(i, j)));
                }

                for (int j = 0; j < b.Columns; j++)
                {
                    result.SetValueAtPosition(new MatrixPosition(i, j + a.Columns),
                                              b.GetValueAtPosition(new MatrixPosition(i, j)));
                }
            }

            return result;
        }

        /// <summary>
        /// Unions the specified matrices. The two matrices have to have the same number of
        /// columns.
        /// </summary>
        /// <param name="a">The first matrix.</param>
        /// <param name="b">The second matrix.</param>
        /// <returns>The union matrix.</returns>
        public static Matrix UnionMatricesHorizontal(Matrix a, Matrix b)
        {
            if (a.Columns != b.Columns)
            {
                throw new ArgumentException("The number of columns have to be even in both matricies.");
            }

            Matrix result = new Matrix(a.Rows + b.Rows, a.Columns);

            for (int i = 0; i < a.Rows; i++)
            {
                for (int j = 0; j < a.Columns; j++)
                {
                    result.MatrixData[i, j] = a.MatrixData[i, j];
                }
            }

            int offset = a.Rows;

            for (int i = 0; i < b.Rows; i++)
            {
                for (int j = 0; j < b.Columns; j++)
                {
                    result.MatrixData[i + offset, j] = b.MatrixData[i, j];
                }
            }

            return result;
        }

        /// <summary>
        /// Implements the operator *.
        /// </summary>
        /// <param name="matrix">The first matrix..</param>
        /// <param name="multiply">The second matrix.</param>
        /// <returns>The result of the operator.</returns>
        public static Matrix operator *(Matrix matrix, Matrix multiply)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (multiply == (Matrix) null)
            {
                throw new ArgumentNullException("multiply");
            }

            if (matrix.Columns == multiply.Rows)
            {
                Matrix result = new Matrix(matrix.Rows, multiply.Columns);
                double[,] resultData = result.MatrixData;
                double[] bcolj = new double[matrix.Columns];

                for (int j = 0; j < multiply.Columns; j++)
                {
                    for (int k = 0; k < matrix.Columns; k++)
                    {
                        bcolj[k] = multiply.MatrixData[k, j];
                    }

                    for (int i = 0; i < matrix.Rows; i++)
                    {
                        double[] arowi = new double[matrix.Columns];

                        for (int l = 0; l < matrix.Columns; l++)
                        {
                            arowi[l] = matrix.MatrixData[i, l];
                        }
                        double s = 0;

                        for (int k = 0; k < matrix.Columns; k++)
                        {
                            s += arowi[k] * bcolj[k];
                        }

                        resultData[i, j] = s;
                    }
                }

                return result;
            }

            throw new IllegalArithmeticException(
                "The number of columns of the one side operator and the number of rows of the other side operator have to be even");
        }

        /// <summary>
        /// Implements the operator +.
        /// </summary>
        /// <param name="matrix">The first matrix.</param>
        /// <param name="add">The second matrix.</param>
        /// <returns>The result of the operator.</returns>
        public static Matrix operator +(Matrix matrix, Matrix add)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (add == (Matrix) null)
            {
                throw new ArgumentNullException("add");
            }

            if ((matrix.Rows == add.Rows) && (matrix.Columns == add.Columns))
            {
                Matrix result = matrix;

                for (int i = 0; i < result.Rows; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        result.SetValueAtPosition(new MatrixPosition(i, j),
                                                  result.GetValueAtPosition(new MatrixPosition(i, j))
                                                  + add.GetValueAtPosition(new MatrixPosition(i, j)));
                    }
                }

                return result;
            }

            throw new IllegalArithmeticException("The matrix dimension has to be even");
        }

        /// <summary>
        /// Implements the operator -.
        /// </summary>
        /// <param name="matrix">The first matrix.</param>
        /// <param name="sub">The second matrix.</param>
        /// <returns>The result of the operator.</returns>
        public static Matrix operator -(Matrix matrix, Matrix sub)
        {
            if (matrix == (Matrix) null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (sub == (Matrix) null)
            {
                throw new ArgumentNullException("sub");
            }

            if ((matrix.Rows == sub.Rows) && (matrix.Columns == sub.Columns))
            {
                Matrix result = matrix;

                for (int i = 0; i < result.Rows; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        result.SetValueAtPosition(new MatrixPosition(i, j),
                                                  result.GetValueAtPosition(new MatrixPosition(i, j))
                                                  - sub.GetValueAtPosition(new MatrixPosition(i, j)));
                    }
                }

                return result;
            }

            throw new IllegalArithmeticException("The matrix dimension has to be even");
        }

        /// <summary>
        /// Implements the operator ==.
        /// </summary>
        /// <param name="matrix">The matrix.</param>
        /// <param name="compare">The compare matrix.</param>
        /// <returns>The result of the operator.</returns>
        public static bool operator ==(Matrix matrix, Matrix compare)
        {
            if (System.Object.ReferenceEquals(matrix, compare))
            {
                return true;
            }

            if (((object) matrix == null) || ((object) compare == null))
            {
                return false;
            }

            return matrix.CompareTo(compare);
        }

        /// <summary>
        /// Implements the operator !=.
        /// </summary>
        /// <param name="matrix">The matrix.</param>
        /// <param name="compare">The compare matrix.</param>
        /// <returns>The result of the operator.</returns>
        public static bool operator !=(Matrix matrix, Matrix compare)
        {
            return !(matrix == compare);
        }

        /// <summary>
        /// Determines whether the current matrix is a square matrix.
        /// </summary>
        /// <value>
        /// True if the current matrix is a square matrix otherwise, false. 
        /// </value>
        public bool IsSquareMatrix
        {
            get { return this.rows == this.columns; }
        }

        /// <summary>
        /// Determines whether the current matrix is symmetric.
        /// </summary>
        /// <value>
        /// True if the current matrix is symmetric otherwise, false.
        /// </value>
        public bool IsSymmetricMatrix
        {
            get
            {
                if (!this.IsSquareMatrix)
                {
                    return false;
                }

                for (int i = 0; i < this.rows; i++)
                {
                    for (int j = 0; j < i; j++)
                    {
                        if (this.MatrixData[i, j] != this.MatrixData[j, i])
                        {
                            return false;
                        }
                    }
                }

                return true;
            }
        }

        /// <summary>
        /// Determines whether the current matrix is a diagonal matrix.
        /// </summary>
        /// <value>
        /// True if the current matrix is a diagonal matrix otherwise, false.
        /// </value>
        public bool IsDiagonalMatrix
        {
            get
            {
                if (!this.IsSquareMatrix)
                {
                    return false;
                }

                for (int i = 0; i < this.rows; i++)
                {
                    for (int j = 0; j < i; j++)
                    {
                        if (i != j)
                        {
                            if (this.matrixData[i, j] != 0)
                            {
                                return false;
                            }
                        }
                    }
                }

                return true;
            }
        }

        /// <summary>
        /// Determines whether the current matrix is a column vector.
        /// </summary>
        /// <value>
        /// True if the current matrix is a column vector otherwise, false.
        /// </value>
        public bool IsColumnVector
        {
            get { return this.columns == 1; }
        }

        /// <summary>
        /// Determines whether the current matrix is a row vector.
        /// </summary>
        /// <value>
        /// True if this instance is row vector otherwise, false.
        /// </value>
        public bool IsRowVector
        {
            get { return this.rows == 1; }
        }

        /// <summary>
        /// Gets the trace of the matrix.
        /// </summary>
        /// <value>The trace of the matrix.</value>
        public double Trace
        {
            get
            {
                double trace = 0;

                for (int i = 0; i < Math.Min(this.rows, this.columns); i++)
                {
                    trace += this.matrixData[i, i];
                }

                return trace;
            }
        }

        /// <summary>
        /// Copies the current instance.
        /// </summary>
        /// <returns>A copy of the current instance.</returns>
        public Matrix Copy()
        {
            return new Matrix((double[,]) this.matrixData.Clone());
        }

        /// <summary>
        /// Checks if the current instance of the matrix is equals
        /// to another.
        /// </summary>
        /// <param name="matrix">The matrix to compare.</param>
        /// <returns>Returns true if the matrix is even to the specified compare matrix else the method returns false.</returns>
        public bool CompareTo(Matrix matrix)
        {
            if ((this.columns != matrix.Columns) || (this.rows != matrix.Rows))
            {
                return false;
            }

            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    if (this.GetValueAtPosition(new MatrixPosition(i, j)) !=
                        matrix.GetValueAtPosition(new MatrixPosition(i, j)))
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        /// <summary>
        /// Determines whether the specified <see cref="T:System.Object"/> is equal to the current <see cref="T:System.Object"/>.
        /// </summary>
        /// <param name="obj">The <see cref="T:System.Object"/> to compare with the current <see cref="T:System.Object"/>.</param>
        /// <returns>
        /// true if the specified <see cref="T:System.Object"/> is equal to the current <see cref="T:System.Object"/>; otherwise, false.
        /// </returns>
        /// <exception cref="T:System.NullReferenceException">The <paramref name="obj"/> parameter is null.</exception>
        public override bool Equals(object obj)
        {
            return base.Equals(obj);
        }

        /// <summary>
        /// Serves as a hash function for a particular type.
        /// </summary>
        /// <returns>
        /// A hash code for the current <see cref="T:System.Object"/>.
        /// </returns>
        public override int GetHashCode()
        {
            return base.GetHashCode();
        }

        /// <summary>
        /// Returns a <see cref="T:System.String"/> that represents the current <see cref="T:System.Object"/>.
        /// </summary>
        /// <returns>
        /// A <see cref="T:System.String"/> that represents the current <see cref="T:System.Object"/>.
        /// </returns>
        public override String ToString()
        {
            String result = "";

            for (int i = 0; i < this.rows; i++)
            {
                for (int j = 0; j < this.columns; j++)
                {
                    result += this.matrixData[i, j] + "\t";
                }

                result += "\n";
            }

            return result;
        }
    }
}